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Abstract 

A method to calculate the crystal field parameters ab initio is proposed and applied to trivalent rare earth 
' impurities in yttrium aluminate and to Tb 3+ ion in TbA103. To determine crystal field parameters local 

O 

^S) ■ Hamiltonian expressed in basis of Wannier functions is expanded in a series of spherical tensor operators. 
5— i ' 

■ Wannier functions are obtained by transforming the Bloch functions calculated using the density functional 

5: 

theory based program. The results show that the crystal field is continuously decreasing as the number of 
4/ electrons increases and that the hybridization of 4/ states with the states of oxygen ligands is important. 

D . Theory is confronted with experiment for Nd 3+ and Er 3+ ions in YAIO3 and for Tb 3+ ion in TbA103 and 

+r| , a fair agreement is found. 

a , 

g . PACS numbers: 71.70.Ch,78.20.Bh,71.15.Mb 
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I. INTRODUCTION 



Ab initio calculations of the properties of molecules and solids became a common tool of solid states 
physics and quantum chemistry. Nevertheless there remain problems which are still not satisfactorily 
treated by the ab initio methods, one of these problems being description of the 4/ states of rare earth 
(R) atoms and ions. Providing that 4/ electrons are localized, the effective Hamiltonian method can 
be applied when analyzing the experimental results. Effective Hamiltonian depends on parameters, 
values of which must be fixed, however. To this end either the fit to experimental results may be 
used or the parameters must be estimated using semiempirical or ab initio methods. 

Original motivation of this work was to explain the physics of magnetic rare earth cobaltites 
RC0O3 (R = rare earth). In these compounds only few experimental data are available, their number 
is certainly insufficient to estimate the crystal field parameters (CFP), hence a necessity to calculate 
CFP emerged and to this end we developed a method described in this paper. To check reliability 
of the method it is applied to orthorhombic rare earth aluminates, which possess the same crystal 
structure as RC0O3. The systems R^Yi^AlOa are widely used for lasers, scintillator and optical 
recording media, rich experimental data exist for them and in several cases reliable set of CFP was 
deduced by fitting the theory and experiment (see e.g.->2). 

Routinely the effective Hamiltonian is written as a sum of the free ion (atomic) part Ha and part 
Hqf describing the crystal field: 

H = H A + H CF . (1) 

The free ion Hamiltonian is spherically symmetrical and it depends on many parameters. Approx- 
imate values of these parameters are either known, or they may be calculated, however. Details 
concerning the operators and parameters of Ha are described e.g. in Ref.-. 

The crystal field Hamiltonian represents more formidable problem. Within a single electron, 
crystal field theory it may be written as^ 

k=0 q=—k i 

where Cg k \i) is a spherical tensor operator of rank k acting on ith electron and the summation 
involving i is over the / electrons of R ion. B q k ^ are crystal field parameters. For the / electrons 
maximum value k max of k is equal to six, providing that cross terms of Hqf between different electron 
configurations are neglected. Hermiticity of Hcf requires B_ = (—l) q B^*. 

(k) 

The values of q and k for which B q are nonzero depend on the site symmetry. If the local 
symmetry of R site is low the number of CFP may be large and it is not possible a priori predict, 
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which of them will be important. When analyzing the experimental results the CFP are therefore 
usually determined from the best agreement of calculated and measured results using the least squares 
method. This is difficult, often ambiguous procedure, and it requires starting values for CFP, however. 
In magnetic or superconducting compounds the number of experimental data is usually insufficient 
to determine CFP, yet magnetic properties reflect the crystal field sensitively. There were therefore 
numerous attempts to estimate CFP theoretically and the effort is continuing (see Ref.- for recent 
survey) . 

In the present paper we propose to use a local Hamiltonian, expressed in the basis of the Wannier 
functions to calculate ab initio the CFP. The method is applied to determine CFP of R 3+ ions 
in orthorhombic aluminates. Comparison with experiment, presented in section |V] shows that the 
method is capable of reliable prediction of the 4/(R) levels splitting by the crystal field. ' 



II. CRYSTAL FIELD STATES OF R 3+ IONS IN ORTHORHOMBIC ALUMINATES 

Rare earth aluminates in question have orthorhombically distorted perovskite structure, belonging 
to the symmetry space group. The unit cell of RAIO3 contains four formula units. R ions are 
located on sites with C s point symmetry and are surrounded by twelve oxygen ligands. 

For R in orthorhombic aluminates a horizontal symmetry plane causes all with odd q to be 
equal to zero, leaving three real B^ parameters (k= 2, 4, 6) and six independent complex parameters 
Bq k ^ (k= 2, 4, 6; q = 2, 4, 6; q ^ k). There are thus fifteen quantities to be determined. It is known& 
that for site symmetries allowing complex CFP, an arbitrary rotation around the orthorhombic c 
axis leaves the eigenvalues of Hamiltonian (P) unchanged. This allows elimination of the imaginary 
part of one of the CFP, easing thus the least squares determination of B^ from the experimental 

results, at the expense that the angle of rotation remains undetermined. In a standard convention 

(2) 

.E?2 parameter is selected to be real. In the method we propose all fifteen quantities are determined. 
The rotation is only invoked in section IVTl when comparing calculated CFP with those obtained using 
least squares by Duan et al— and Gruber et al.—. 

Seven 4/ states \l,m) (1=3, m = ±3, ±2, ±1,0) of the R are split in C s symmetry crystal field 
in seven singlets. The singlets are of two different types. Taking the axis of quantization along the 
orthorhombic c axis four of these singlets are formed by m = ±3 and m = ±1 states, in remaining 
three singlets ±2 orbitals are mixed with m =0 state. In our analysis instead of the \l,m) states the 
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basis of real orbitals is used 
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In terms of real orbitals the wave functions of the seven singlets may be written as 

4 7 
tpi = ^2cj,i\<Pj); z = l,2,3,4; il>i = ^2cj,i\<Pj); i = 5,6,7. (4) 

j=l j=5 

III. CALCULATION OF ELECTRONIC STRUCTURE AND DESCRIPTION OF METHOD 

Starting step of our analysis is a standard calculation of the electronic structure of compounds 
in question. To this end we used the augmented plane waves + local orbital method based on the 
density functional theory as implemented in the WIEN2k program^. For the exchange-correlation 
functional the generalized-gradient approximation form& was adopted. 

For the orthorhombic unit cell parameters their experimental values of YAIO^ and TbAlO^ were 
taken, but the internal crystal structure parameters were optimized by minimizing the forces acting 
on the atoms. Typical content x of R ions in Yi^R-cAlC^ used as laser material is x ~ 0.01-0.03. In 
our calculations the unit cell contained 120 atoms (RY23AI24O72, corresponding to x=0.0435) and it 
retains orthorhombic symmetry. The eigenvalue problem was solved in four points of the irreducible 
Brillouin zone and number of basis functions was ~ 9200 (corresponding to parameter RK max =6. 13). 
The calculations were nonspinpolarized and the 4/ electrons of R were treated as core electrons. 

Core states contribute to the spherically symmetrical component of the density only and, as a 
consequence, the potential on the lattice site of R in question does not contain any nonspherical 
component arising from the on-site R 4/ states. This is vital when determining the crystal field, as 
otherwise nonphysical interaction of the 4/ states with nonspherical potential they create themselves 
(selfinteraction) would give dominating contribution to CFP. 
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In the second step the 4/(R) states are included in the set of valence states and the eigenvalue 
problem is solved with the potential obtained in calculation described above. The goal is to obtain 
reliable density of states projected on the 4/ subspace and to determine from the 4/ Bloch states 
corresponding Wannier functions. Difficult problem is the hybridization of 4/ states with other 
valence states. It depends on relative position of the energy levels and no available approximation 
to DFT gives the position of 4/ levels correctly. To deal with this problem we shift the levels using 
on-site orbital potentials. For atom i and state a (a = s,p, d, f) the form of such potential is 

where A$ )Q is the magnitude of the shift and P^ a is a projection operator on states i, a. In aluminates 
the crucial is position of 4/ levels relative to the oxygen 2p and to a lesser extent 2s levels and this 
makes the results of our calculation dependent on a parameter 



*ox,2p,2s- 



(6) 



The procedure indicating how A can be estimated is described in the Appendix. We checked that the 
results are only slightly influenced by a reasonable shift of other levels (see section HV|) . To remove 
any ambiguity, the energy of oxygen levels was shifted by A, no shift was applied to 4/(R) states, 
while energy of all remaining valence states was increased by 6 Ry (~ 82 eV) preventing thus their 
hybridization with 4/(R) orbitals. 

The distance between the nearest R ions in RY23AIO72 structure, we used is larger than 1 nm, 
hence the hybridization of 4/ states on different lattice sites is small and it is not difficult to establish 
correspondence between the peaks in 4/ DOS and the 4/ energy levels. In TbAlC>3 the situation is 
more complex. The Tb-Tb distance is only 0.365 nm and 4/-4/ hybridization leads to a considerable 
dispersion of the 4/ bands, which makes the analysis of the DOS complicated. To circumvent the 
problem we made the four Tb in the unit cell of TbA103 formally inequivalent and left 4/ states of 
three of them in the core. Even then there exists a residual dispersion, which makes the DOS peaks 
considerably broader compared to the DOS of Yx^ThrAlOs. The situation is documented in Fig. 

m 

The density of states projected on 4/ states makes possible to determine the energies e$ of 4/ 
eigenstates and also their orbital composition i.e. modules \cij\ of expansion coefficients in eq. fllj. 
£j and \cij\ may be then used to determine the crystal field on the R site. Note that if the peaks in 
the density of states overlap, a decomposition of the peaks on their components must precede the 
analysis. Such a strategy was successfully used to expose the importance of hybridization for Pr 4+ 
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FIG. 1. (Color online) TbAlC-3 and Tb:YA103. Density of states projected on the 4/ subspace. The shift 
A equals to -8.2 eV. E c j is the energy difference between the highest and lowest 4/ singlet state and it is 
used to characterize the strength of the crystal field. 

ion on sites with the cubic symmetry in PrC"2 and PrBaO^. For cubic local symmetry the energies 
Si are sufficient to determine unambiguously the CFP, but if the symmetry is lower, the \cij\ must 
be employed as well. Our attempt to determine CFP in orthorhombic aluminates using \cij\ and 
the least squares method led to several sets of results, however. The situation is reminiscent of the 
determination of CFP from optical spectra (see e.g. Ref.— ). 

In this paper we propose an elegant and unambiguous way to determine the CFP employing the 
Wannier functions. Wannier functions are real space counterpart of the Bloch functions and both sets 
are connected by a unitary transformation. In order that this transformation is unique an additional 
requirement on the form of Wannier functions is needed. One such requirement - maximal localization 
is implemented in the wannier90 package^. The interface between the WIEN2k and wannier90 was 
recently written by Kunes et al.— . Among other outputs wannier90 yields also Hamiltonian H\y in 
basis of the Wannier functions 

H w = U + e k U, (7) 

where U is the matrix of transformation and e& are energies of the Bloch eigenf unctions. From H\y 
we extract the local Hamiltonian H^f - the seven by seven matrix consisting of the matrix elements 
between the 4/ Wannier functions centered on the R ion in question. The local Hamiltonian has 
the symmetry of the R site and it can be expressed as a combination of the unit operator I and the 
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crystal field Hamiltonian fl2]) 

H 4f = E avg I + H CF = E avg I + B™CW(i). (8) 

k,q,i 

In the subspace of the 4/ Wannier functions 7 and Cq (i) may be represented by 7x7 matrices and 

(k) 

the above equality yields the system of linear equations from which E avg and B q are determined. 

The method was first tested by determining the crystal field parameters of Pr 4+ ion in PrO^. The 
results were similar to those reported by us earlier—, albeit without the uncertainty in decomposition 
of the density of states peaks. We note that concept of the local Hamiltonian is also in the heart of 
CFP calculation method proposed recently by Hu et al.—. 

IV. RESULTS 

The calculation with the 4/(R) states included between the core states resulted in the band 
structure corresponding to an insulator with the band gap ~ 6 eV. The valence band is dominated 
by oxygen 2p states, while the bottom of the conduction band is mainly formed by 5(i(R) and 4d(Y) 
orbitals. 

When the 4/(R) states are treated as valence states (second step in our method) the results depend 
naturally on the specific R 3+ ion and on the local geometry, but also on the value of parameter A. As 
|A| decreases hybridization of 4/(R) with oxygen 2p states becomes larger and for |A| smaller than 
~ 2.7 eV the gap between these levels disappears and the 4/ density of states becomes unresolved. 
The strength of the crystal field may be characterized by the difference E c f of the lowest and highest 
4/ eigenenergy (see Fig. 1). E c f for all thirteen R 3+ ions in YAIO3 and for several values of A is 
shown in Fig. 2. As seen from this Figure the dependence of E c f on the number N^f of 4/ electrons 
is smooth and monotonically decreasing. The hybridization represented by the shift A plays an 
important role. 

The nonzero CFP of R 3+ ions in YAIO3 are displayed in Figs. [3] and H] for typical value of the 
shift A= -8.2 eV. The dependence of B k q on JV4/ is again smooth, the largest energy is carried by the 
k=6, q = 4 term. 

Returning to Fig. [TJ we note that the crystal field in TbAlOs is markedly stronger compared to 
the one in Tb:YA103. This indicates a large effect of the local geometry - in both cases this geometry 
was determined by minimizing the forces on atoms constituting the system in question. 

We remind that when obtaining above described results, the 4/(R) states were allowed to hybridize 
only with oxygen 2p and 2s states - hybridization with other valence states was prevented by shifting 
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FIG. 2. (Color online). R:YA103. Difference E c f of the lowest and the highest 4/ eigenenergy on number 
of the 4/ electrons for several values of the shift A. Circles correspond to a calculation in which all valence 
states, except 4/(R) were shifted by 81.6 eV (6 Ry), preventing any hybridization of the 4/ states. The 
curves in this, as well as in the following fieures. serve as euides for eves onlv. 
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FIG. 3. (Color online). R:YA103. Dependence of real part of crystal field parameters on number of the 4/ 
electrons. Shift A= -8.2 eV. 
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FIG. 4. (Color online). R:YA103. Dependence of imaginary part of crystal field parameters on number of 
the 4/ electrons. Shift A= -8.2 eV. 

their energy by ~82 eV. For Er:YA103 system we analyzed in details how the results are influenced by 
such an approximation. The eigenvalue problem was solved with A=-8.2 eV and energy of selected 
type of valence state left unshifted, allowing thus its hybridization with 4/(R). The importance of 
the hybridization may be characterized by quantity Ae^ 

ASi = ( £i - £?)/£?, (9) 

where is energy of ith 4/ level for all Er, Y and Al states shifted high in energy, while E{ corresponds 
to calculation in which selected valence state was not shifted. The results collected in Table [I] show 
that hybridization of 4/ states with other valence states of the cations changes the 4/ energies by 
few %. The most pronounced effect has the hybridization with 5c? states of erbium. 



V. COMPARISON WITH EXPERIMENT 

In Fig. [5] the energies of Tb 3+ seven lowest multiplets are displayed. The energies are taken 
relative to the lowest energy of the multiplet in question. Experimental data£ are compared with 
results calculated for two values of A, which give the smallest deviation from the experiment. 

A similar plot for the Nd 3+ in YAIO3 is presented in Fig. |6j In this case the results are more 
sensitive to the value of A and calculations for A = -4.1, -8.2 and -10.9 eV are confronted with the 
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FIG. 5. (Color online)TbA103. Energies of Tb 3+ eigenstates taken relative to the lowest energy of the 
multiplet. The experimental data were determined by Gruber et al£. 

experiments^. In addition the data obtained with all, but 4/ states displaced by 6 Ry (~ 82 eV), 
which prevented the hybridization, are also included. 

As the third example we consider Er 3+ ion in YAIO3. As seen in Fig. [7] the experimental data 
for the four lowest multiplets are in very good agreement with the calculation in which the shift 
parameter A is equal to -8.2 eV. For Er:YAlC>3 and A = -8.2 eV we calculated the mean square 
deviation a of the experimental and calculated splittings: 

"'exp j—-y 
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n (index of eigenstate) 



FIG. 6. (Color online)Nd:YA103. Energies of Nd 3+ eigenstates taken relative to the lowest energy of the 
multiplet. The experimental data were determined by Kaminskii2£ and Duan et al.—. 

where n exp is the number of experimentally observed levels in the jth set and the energies are again 
taken relative to the lowest level in the set (the sets correspond to the \L, S, J, Mj) multiplets, with 
the exception of the set 14, which combines 2 .^13/2, 2 Pi/2 and A G§/2 multiplets). In Fig. |S]the mean 
square deviation is displayed for the lowest 19 Er 3+ sets. With the exception of sets 14, 16 and 19 
values of a indicate a fair agreement of experiment and calculation. We note that in these three 
cases not all levels were observed experimentally, while for the remaining 16 multiplets all levels were 
detected. 



VI. DISCUSSION 

The results presented in Figs. 2, 6 and 7 show convincingly that the hybridization of 4/ states is 
important. Our semiempirical estimation of the parameter A, which reflects its strength, amounts 
to -4.3, -11.5 and -5.1 eV for Nd, Tb and Er, respectively (see the Appendix). This qualitatively 
agrees with A for which the calculated positions of the energy levels fit the experiment best (optimal 
|A| being in limits 5-8, 8-11 eV for Nd and Tb and ~ 8 eV for Er). 

In order to compare calculated CFP with those obtained by least squares fit to the optical 
spectral, the rotation angle a must be first calculated, as explained in section [Til It is straight- 
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FIG. 7. (Color online)Er:YA103. Energies of Er 3+ eigenstates taken relative to the lowest energy of the 
multiplet. Thi — 
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FIG. 8. (Color online)Er:YA103. The mean square deviation a of the experimental splittings and splittings 
calculated with A= -8.2 eV. Crosses correspond to the multiplicity of the set of levels, number of levels 
observed experimentally n exp is denoted by squares. The experimental data were determined by Donlan 
and Santiago^ and Duan et al.—. 
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FIG. 9. (Color online)Angle of rotation around orthorhombic c axis. 

forward to show that if imaginary part of B\ in rotated system is to be zero, a is given by 

a = -arctg(Bl imag /Bl real ) , (11) 

where B\ imag , B\ Teal are imaginary and real parts of k—2, q=2 CFP in orthorhombic coordinate 
system. As seen from Figs. El H] for A = -8.2 eV B\ imag is large and positive for all R 3+ ions, 
while B\ real is small and positive for lighter R, changing the sign for R = Dy. The dependence 
of a on N^f is displayed in Fig. |9j As seen from this Figure a(N±f) change the sign, which leads 
to a discontinuity of B^(N^f) dependence despite the fact that in original orthorhombic system the 
crystal field parameters change continuously with N^f. This concerns all q=2 and q =4 terms and 
for this reason it is more informative to compare absolute values of B^, rather than to compare 
separately their real and imaginary parts. Such comparison for Nd 3+ and Er 3+ ions in YAIO3 and 
for Tb 3+ ion in TbA103 is given in Table [TTJ The correspondence between calculated CFP and those 
obtained by the least squares fit of the experiment is satisfactory. In particular the k = 6, q = 4 
term carries the largest energy and the crystal field decreases with increasing number of 4/ electrons. 
CFP determined by fitting the optical spectra suffer from the ambiguity connected with numerous 
local minima of the least squares method applied in the 14 dimensional space of variables. In this 
connection we return to the Er:YA103 case (Fig. |8]) in which the agreement between our calculation 
and the energy levels fixed using optical spectra was significantly worse for the level sets 14, 16, 19, 
for which not all the levels were detected. The assignment of measured levels in the least squares 
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TABLE II. Comparison of absolute values of crystal field parameters calculated with the shift A=-8.2 eV 
with those obtained by Duan et a/ji (a) and Gruber et al£ (b) by least squares fit to optical spectra. All 
CFP are in units of cm . 

method is far from being unique and the results of our calculation may serve as a useful seed values 
for it. 

There are several limitations of the method we suggest for the CFP calculation. In order to prevent 
the selfinteraction, the 4/ hybridization with other valence states was treated in a nonselfconsistent 
way. This may be a problem for Ce 3+ ion and for compounds with tetravalent R ions, where the 
hybridization is known to be significant. For trivalent R heavier than Ce, the effect is likely to be 
small. 

Another approximation is connected with the assumed single electron character of the crystal 
field Hamiltonian. Obviously the closer \L, S, J, Mj) multiplet is to the top of the valence band, the 
stronger is the hybridization. This leads to CFP being different in different multiplets, the effect 
which is included if the concept of correlation crystal field^ 3 ^ is used, or if CFP are allowed to be 
energy dependent. Taking into account that different multiplets of Nd 3+ , Tb 3+ and Er 3+ are well 
described by single CFP set (section EJ), CFP energy dependence and correlation crystal field play 
only a minor role in cases considered here. 

While the above discussed crystal field Hamiltonian yields satisfactorily the positions of 4/ levels, 
it gives no information on the intensity of 4/ — 4/ transitions and the Judd-Ofelt theory^ must be 
used to this end. It was suggested recently by Hu et alM that concept of 'local dipolar operator', 
which is similar to the concept of local Hamiltonian, may be used instead of the Judd-Ofelt theory. 
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FIG. 10. (Color online)Scheme of levels in rare earth containing oxides. 

Very recently a part, which allows to treat optics, was added to the wien2wannier package^ 3 - 1 ^. It 
is relatively straightforward to extract from that part the 'local dipolar operator' expressed in the 
basis of Wannier functions. The generalization of the method we propose here will then allow to 
determine both positions of 4/ levels and 4/ — 4/ transition intensities in a single step. 

VII. CONCLUSIONS 

In conclusion we presented a method to calculate the CFP that is fully ab initio except a single 
parameter A which corresponds to energy difference between the unsplit 4/(R) and oxygen valence 
levels. A is estimated using independent experimental data and its estimation may be useful when 
explaining other experiments e.g. the luminescence. The viability of the method is demonstrated by 
the good agreement between the calculated and experimentally observed crystal field splitting of the 
\L, S, J, Mj) multiplets of Nd 3+ , Tb 3+ and Er 3+ ions. The method may become a valuable tool to 
interpret the experiments, especially after the determination of 4/ — 4/ transition probabilities will 
be added as described in section IVT1 

Appendix A: Energy of rare earth 4/ level relative to oxygen 2p states 

The scheme of energy levels in the rare earth aluminates is shown in Fig. [TOj As seen from this 
Figure the shift A may be written as a sum of three contributions 

A = 5 L + S D + 6 G , (Al) 
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FIG. 11. (Color online) Total density of states of TDY23AI24O72 with 4/ states included between the core 
states. Dashed line corresponds to the center of gravity of the valence band, Fermi energy is at zero. 

where 5l is difference of unsplit 4/ level energy and the energy of the ground \L, S, J, Mj) multiplet. 
8d is difference of the ground multiplet energy and energy of the top of valence band. Finally 8q is 
obtained by subtracting energies of the top and center of gravity of this band. 

The valence band is dominated by the oxygen 2p states, width of this band in the orthorhombic 
aluminates is ~ 6.9 eV and the energy difference 5g is ~ 1.84 eV (Fig. [IT]) . 

In a series of papers Dorenbos collected and analyzed the spectroscopic data concerning the 
4/ electron states of rare earth substitutions in transition metal complex compounds (see^ and 
references therein). In particular the energy difference, we denote as 5d was determined. Sd depends 
on the compound in question, but as discussed by Rogers et al.— the difference #d(R) - #d(R') is to a 
good approximation compound independent. In Fig. H2]^z)(R) - 5o(Gd) is displayed for all trivalent 
rare earth ions. The data were adopted from Fig. 1 of reference^. In seven different oxides for which 
Dorenbos et al.— analyzed the spectroscopic data the ground 8 SV/2 multiplet of Gd 3+ lies below the 
top of the valence band, 5o(Gd) is thus negative with values ranging from -5 to -2.5 eV. 

To determine the energy difference Sl between the unsplit (4/) n and the lowest \L, S, J, Mj) 
multiplet we used the 'lanthanide' program^ with parameters of Carnall et al.— . The dependence 
of 5l on the type of trivalent rare earth is shown in Fig. fT2l 

Finally adding <5d(R) - <5o(Gd), 5l(R), 5q with minimal and maximal values of 5o(Gd) yields the 
range in which A(R) should lie (Fig. IT3|) . 
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FIG. 12. (Color online) Energy differences <5l(R) between the unsplit 4/ level and the ground multiplet 
(open circles), and <$d(R) of the ground multiplet and the top of the valence band. The data for <5d(R) are 
relative to ^(CnM and were adonted from refji^. 
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FIG. 13. (Color online) Minimum and maximum of the energy difference A between the unsplit 4/ level 
and the center of gravity of the valence band, n is number of electrons of trivalent rare earth ion. 
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